SUBROUTINE INMSOIL_INIT &
& (tsnoi,wgnoi,xdd2,snowini,dt)

use MPI_VARIABLES, only : &
& nxi, nx1e, nyi, ny1e, &
& nxs2i, nxs2e, nys2i, nys2e

use INMSOIL_ARRAYS

use INMSOIL_INMSOIL

use INMSOIL_NUMPARAMS, only : &
& KL, ML, MS, NT, Num_Soil, &
& Num_Veget, dznorm, dzmin, UpperLayer

use INMSOIL_PHYSPARAMS, only : &
& POR_soil

use ALLOC, only : &
& Kelvin0

implicit none

include 'soil.inc'

! Input variables

real(8), intent(in) :: tsnoi(nxs2i:nxs2e,nys2i:nys2e) ! Initial soil
                                                      ! temperature, Kelvin
real(8), intent(in) :: wgnoi(nxs2i:nxs2e,nys2i:nys2e) ! Initial soil moisture,
                                                      ! kg/kg
real(8), intent(in) :: xdd2(nxs2i:nxs2e,nys2i:nys2e)  ! The depth of soil layer, m
                                                      ! (domain)
real(8), intent(in) :: snowini(nxs2i:nxs2e,nys2i:nys2e)  ! Initial depth of snow, m
                                                      
real(8), intent(in) :: dt

! Local variables

REAL, allocatable :: T_INI(:) ! ML
REAL :: Q_INI, GSOLO_INI
REAL :: WSOIL_INI, WSKIN_INI
REAL, allocatable :: WSG_INI(:), VSG_INI(:), WIG_INI(:) ! ML
REAL, allocatable :: Q_METHANE_INI(:) ! ML
REAL :: PRECIP_INI
REAL :: SN_INI, SN_INTRCP_INI 
REAL, allocatable :: SN_DZ_INI(:), SN_DENSITY_INI(:) ! MS

real(8), save :: pres0 = 1.d+5 ! Typical surface atmospheric pressure, Pa
real(8), save :: sndens0 = 0.2 ! Typical snow density, g/cm**3i
real(8), parameter :: m_to_cm = 1.d+2

integer(4) :: ix, iy ! Loop indices
integer(4) :: k

real(8), external :: QSAT

!logical, save :: firstcall = .true.

call INMSOIL_ARRAYS_ALLOC(nxi,nx1e,nyi,ny1e,ml,ms)

! assigning boundary layer parameterization constants
call PBLDAT_SOIL( &
& FBL,vkc,anu,z0min,chc,zeta_st, &
& p_m,p_t,q_m,q_t,alpha_m,alpha_t, &
& beta_m,beta_t,gamma_m,gamma_t, &
& a_m,b_m,a_t,b_t )

allocate (T_INI(ML))
allocate (WSG_INI(ML), VSG_INI(ML), WIG_INI(ML))
allocate (Q_METHANE_INI(ML))
allocate (SN_DZ_INI(MS), SN_DENSITY_INI(MS))

do iy = nyi, ny1e
  do ix = nxi, nx1e
    PRECIP_ACC(ix,iy)=0.0
    T_INI(1:ml) = tsnoi(ix,iy) ! Kelvin degrees
    Q_INI = QSAT(tsnoi(ix,iy),pres0) ! kg/kg
    GSOLO_INI = 0.
!    WSOIL_INI = wgnoi(ix,iy)*xdd2(ix,iy)*100. ! 100 - transform to centimeters
    WSOIL_INI = wgnoi(ix,iy)*(1.-POR_soil)*xdd2(ix,iy)*1.2*m_to_cm ! 1.2 - ratio of soil to water density;
                                                                   ! for soiltype = 0 only (!)
    WSKIN_INI = 0.5 !
    do k = 1, ml
      if (T_INI(k)-Kelvin0 > 0.) then
        WSG_INI(k) = wgnoi(ix,iy)
        WIG_INI(k) = 0.
      else
        WSG_INI(k) = 0. !0.5*wgnoi(ix,iy)
        WIG_INI(k) = wgnoi(ix,iy) !0.5*wgnoi(ix,iy)
      endif
    enddo
    VSG_INI(1:ml) = 0. ! zero water vapour content in soil as inital condition
    Q_METHANE_INI(1:ml) = 0. ! initial methane concentration is zero, assuming
                             ! that it will be produced up to saturation in short term after
                             ! starting the integration
    PRECIP_INI = 0. ! no atmospheric precipitation at initial time moment

!   Snow variables group
    SN_DZ_INI(1:ms) = 0.
    SN_DENSITY_INI(1:ms) = 0.
    if (snowini(ix,iy) > 0.) then
      k = int(snowini(ix,iy)*m_to_cm/dznorm)
      SN_DZ_INI(ms-k+1:ms) = dznorm
      SN_DZ_INI(ms-k) = &
      & max(snowini(ix,iy)*m_to_cm - sum(SN_DZ_INI),dzmin)
      SN_DENSITY_INI(ms-k:ms) = sndens0 
    endif            
    SN_INTRCP_INI = 0.
    SN_INI = snowini(ix,iy)*m_to_cm*sndens0/1. ! 1. - water density, g/cm**3

    dc6(1:nv2,ix,iy) = 0.
    cpool_soil(1:nsp,ix,iy) = 0.
    npool_soil(1:nsp,ix,iy) = 0.

    call INITSOIL( &
    & T_INI, &         ! snow and soil collective initial temperature
                  ! (ML profile array)
    & Q_INI, &         ! initial air humidity at the surface (scalar), kg/kg
    & GSOLO_INI, &     ! initial soil flux (0)
    & WSOIL_INI, &     ! integral amount of water in soil (0.75 * 15.6), cm
    & WSKIN_INI, &     ! amount of water to capture by vegetation (0.5)
    & WSG_INI, &       ! amount of water in soil
                  ! (ML profile array)
    & VSG_INI, &       ! amount of water vapor in soil
    & WIG_INI, &       ! amount of ice in soil
                  ! (ML profile array)
    & Q_METHANE_INI, & ! initial methane amount
                  ! (ML profile array)
    & PRECIP_INI, &    ! initial precipitation
    & SN_INI, &        ! initial snow water-equivalent depth
    & SN_INTRCP_INI, & ! snow to capture by pines (0)
    & SN_DZ_INI, &     ! thicknesses of snow layers
    & SN_DENSITY_INI, &! densities of snow layers
   
    & KL, ML, MS, NT, NUM_SOIL, NUM_VEGET, &
    & sngl(DT), DZNORM, DZMIN, sngl(xdd2(ix,iy)), UPPERLAYER, &
     
    & ST(ix,iy),PGR(ix,iy),TGROLD(ix,iy),QGROLD(ix,iy),RADIAT(ix,iy), &
    & PRECIP(ix,iy),WSOLD(ix,iy),SNOLD(ix,iy), &
    & ZS(ix,iy),THSOIL(ix,iy),WHSOIL(ix,iy),BOLD(ix,iy),RF1(ix,iy), &
    & RF2(ix,iy),SNMELT(ix,iy),SNMELT_OLD(ix,iy), &
    & HSNOW(ix,iy),HS(ix,iy),ES(ix,iy),TGRNEW(ix,iy),QGRNEW(ix,iy), &
    & WSNEW(ix,iy),SNNEW(ix,iy),RUNOF1(ix,iy), &
    & ELATOLD(ix,iy),HFOLD(ix,iy),PRS(ix,iy),SHORTWAVE(ix,iy), &
    & BW(ix,iy),BMELT(ix,iy),PRSOLD(ix,iy), &
  
    & WSOIL(ix,iy),WSKIN(ix,iy),T0(ix,iy), &
     
    & TSG(1,ix,iy),WSG(1,ix,iy),VSG(1,ix,iy),WIG(1,ix,iy),BG(ix,iy), &
     
    & Q_METHANE(1,ix,iy),TGRANN(1,ix,iy),FPLANT(ix,iy), &
    & FEBUL(ix,iy),FDIFF(ix,iy),FTOT(ix,iy), &
    & PLANT_SUM(ix,iy),BULL_SUM(ix,iy),OXID_SUM(ix,iy), &
    & RPROD_SUM(ix,iy),ANOX(ix,iy),WAT_TAB(ix,iy), &
     
    & TSRFOL(1,ix,iy),TSRF(1,ix,iy), &
     
    & RHOGR(1,ix,iy),WLM0(1,ix,iy),WLM7(1,ix,iy),BH(1,ix,iy),PSIMAX(1,ix,iy), &
    & POR(1,ix,iy),FLWMAX(1,ix,iy),DLMAX(1,ix,iy),DZ(1,ix,iy), &
     
    & EXTINCT(1,ix,iy),DENSNY(ix,iy),SN_INTRCP(ix,iy), &
     
    & AZ0P(ix,iy),VEG(ix,iy),WL(1,ix,iy),ALBEDO_GR(ix,iy), &
   
    & AL(1,ix,iy),DLT(1,ix,iy),DVT(1,ix,iy),ALL(1,ix,iy),DL(1,ix,iy), &
    & ALV(1,ix,iy),DV(1,ix,iy),Z(1,ix,iy),T(1,ix,iy), &
    & WV(1,ix,iy),WI(1,ix,iy),DENS(1,ix,iy))


  enddo
enddo

deallocate (T_INI)
deallocate (WSG_INI, VSG_INI, WIG_INI)
deallocate (Q_METHANE_INI)
deallocate (SN_DZ_INI, SN_DENSITY_INI)

!if (firstcall) firstcall = .false.
END SUBROUTINE INMSOIL_INIT 
